Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC

Produced in R version 3.2.1 using pomp version 0.72.2.

Introduction: What is Simulation-based Inference?

Overview of our ambitions

Noisy clockwork: Time series analysis of population fluctuations in animals

Six problems of Bjornstad and Grenfell (Science, 2001)

  • Obstacles for ecological modeling and inference via nonlinear mechanistic models:

  1. Combining measurement noise and process noise.

  2. Including covariates in mechanistically plausible ways.

  3. Continuous time models.

  4. Modeling and estimating interactions in coupled systems.

  5. Dealing with unobserved variables.

  6. Modeling spatial-temporal dynamics.

  • The same issues arise for epidemiological modeling and inference via nonlinear mechanistic models.

Partially observed Markov process (POMP) models

  • Data \(y^*_1,\dots,y^*_N\) collected at times \(t_1<\dots<t_N\) are modeled as noisy and incomplete observations of a Markov process \(\{X(t), t\ge t_0\}\).

  • This is a partially observed Markov process (POMP) model, also known as a hidden Markov model or a state space model.

  • \(\{X(t)\}\) is Markov if the history of the process, \(\{X(s), s\le t\}\), is uninformative about the future of the process, \(\{X(s), s\ge t\}\), given the current value of the process, \(X(t)\).

  • If all quantities important for the dynamics of the system are placed in the state, \(X(t)\), then the Markov property holds by construction.

  • POMP models can include all the features desired by Bjørnstad and Grenfell (2001).

Notation for partially observed Markov process models

  • Write \(X_n=X(t_n)\) and \(X_{0:N}=(X_0,\dots,X_N)\). Let \(Y_n\) be a random variable modeling the observation at time \(t_n\).

  • The one-step transition density, \(f_{X_n|X_{n-1}}(x_n|x_{n-1};\theta)\), together with the measurement density, \(f_{Y_n|X_n}(y_n|x_n;\theta)\) and the initial density, \(f_{X_0}(x_0;\theta)\), specify the entire joint density via

\[f_{X_{0:N},Y_{1:N}}(x_{0:N},y_{1:N};\theta) = f_{X_0}(x_0;\theta)\,\prod_{n=1}^N\!f_{X_n | X_{n-1}}(x_n|x_{n-1};\theta)\,f_{Y_n|X_n}(y_n|x_n;\theta).\]

  • The marginal density for sequence of measurements, \(Y_{1:N}\), evaluated at the data, \(y_{1:N}^*\), is

\[ f_{Y_{1:N}}(y^*_{1:N};\theta)=\int f_{X_{0:N},Y_{1:N}}(x_{0:N},y^*_{1:N};\theta)\, dx_{0:N}. \]

To think algorithmically, we define some function calls:

  • rprocess( ): a draw from \(f_{X_n|X_{n-1}}(x_n| x_{n-1};\theta)\)

  • dprocess( ): evaluation of \(f_{X_n|X_{n-1}}(x_n| x_{n-1};\theta)\)

  • rmeasure( ): a draw from \(f_{Y_n|X_n}(y_n| x_n;\theta)\)

  • dmeasure( ): evaluation of \(f_{Y_n|X_n}(y_n| x_n;\theta)\)

What does it mean for methodology to be simulation-based?

  • Oftentimes, simulating random processes is easier than evaluating their transition probabilities.

  • In other words, we may be able to write rprocess() but not dprocess().

  • Simulation-based methods require the user to specify rprocess() but not dprocess().

  • Plug-and-play, likelihood-free and equation-free are alternative terms for simulation-based.

  • Much development of simulation-based statistical methodology has occurred in the past decade.

Practicum: The pomp R package for POMP models

Objectives

  1. To introduce students to the pomp package
  2. To explain how the components of a POMP model are encoded in this package
  3. To give some experience in the use and manipulation of pomp objects

Intro to pomp

  • pomp is an R package for data analysis using partially observed Markov process (POMP) models.

  • Note the distinction: lower case pomp is a software package; upper case POMP is a class of models.

  • pomp builds methodology for POMP models in terms of arbitrary user-specified rprocess(), dprocess(), rmeasure(), and dmeasure() functions.

  • Following modern practice, most methodology in pomp is simulation-based, so does not require specification of dprocess().

  • pomp has facilities to help construct rprocess(), rmeasure(), and dmeasure() functions for model classes of epidemiological interest.

  • pomp provides a forum for development, modification and sharing of models, methodology and data analysis workflows.

The following diagrams show the structure of a POMP model schematically.

state_space_diagram1

Schematic of the structure of a POMP showing causal relations and direction of information flow. The key perspective to keep in mind is that the model is to be viewed as the process that generated the data.

state_space_diagram2

Another POMP model schematic, showing dependence among model variables. State variables, \(x\), at time \(t\) depend only on state variables at the previous timestep. Measurements, \(y\), at time \(t\) depend only on the state at that time.

Example: the Ricker model

Download a script containing the R commands below.

The deterministic skeleton

The Ricker map describes the deterministic dynamics of a simple population: \[N_{t+1} = r\,N_{t}\,\exp(-N_{t})\]

Here, \(N_t\) is the population density at time \(t\) and \(r\) is a fixed value (a parameter), related to the population’s intrinsic capacity to increase. \(N\) is a state variable, \(r\) is a parameter. If we know \(r\) and the initial condition \(N_0\), the deterministic Ricker equation predicts the future population density at all times \(t=1,2,\dots\). We can view the initial condition, \(N_0\) as a special kind of parameter, an initial-value parameter.

Process noise

We can model process noise in this system by making the growth rate \(r\) into a random variable. For example, if we assume that the intrinsic growth rate is log-normally distributed, \(N\) becomes a stochastic process governed by \[N_{t+1} = r\,N_{t}\,\exp(-N_{t}+\varepsilon_{t}), \qquad \varepsilon_{t}\;\sim\;{\mathrm{Normal}\left(0,\sigma\right)},\] where the new parameter \(\sigma\) is the standard deviation of the noise process \(\varepsilon\).

Measurement error

Let’s suppose that the Ricker model is our model for the dynamics of a real population. However, we cannot know the exact population density at any time, but only estimate it through sampling.

Let’s model measurement error by assuming the measurements, \(y_t\), are Poisson with mean \(\phi\,N_t\): \[y_{t}\;\sim\;{\mathrm{Poisson}\left(\phi\,N_{t}\right)}\]

In this equation,

  1. \(N_t\) is the true population density at time \(t\),
  2. \(y_t\) is the number of individuals sampled at time \(t\),
  3. the choice of units for \(N\) is peculiar and depends on the parameters (e.g., \(N=\log(r)\) is the equilibrium),
  4. the parameter \(\phi\) is proportional to our sampling effort, and also has peculiar units.

Working with the Ricker model in pomp.

The R package pomp provides facilities for modeling POMPs, a toolbox of statistical inference methods for analyzing data using POMPs, and a development platform for implmenting new POMP inference methods. The basic data-structure provided by pomp is the object of class pomp, alternatively known as a pomp object. It is a container that holds real or simulated data and a POMP model, possibly together with other information such as model parameters, that may be needed to do things with the model and data.

Let’s see what can be done with a pomp object. First, if we haven’t already, we must install pomp and pompExamples. Run the following to do so.

require(devtools)
install_github("kingaa/pomp")
# install_github("kingaa/pompExamples")

Now we’ll load some packages, including pomp.

set.seed(594709947L)
require(ggplot2)
require(plyr)
require(reshape2)
require(pomp)
stopifnot(packageVersion("pomp")>="0.69-1")

A pre-built pomp object encoding the Ricker model comes included with the package. Load it by

pompExample(ricker)
## newly created object(s):
##  ricker

This has the effect of creating a pomp object named ricker in your workspace. We can plot the data by doing

plot(ricker)

We can simulate by doing

x <- simulate(ricker)

What kind of object have we created?

class(x)
## [1] "pomp"
## attr(,"package")
## [1] "pomp"
plot(x)

Why do we see more time series in the simulated pomp object?

We can turn a pomp object into a data frame:

y <- as.data.frame(ricker)
head(y)
##   time   y
## 1    0  68
## 2    1   2
## 3    2  87
## 4    3   0
## 5    4  12
## 6    5 174
head(simulate(ricker,as.data.frame=TRUE))
##   time  y          N            e sim
## 1    0 68 7.00000000  0.000000000   1
## 2    1  3 0.26614785 -0.069613447   1
## 3    2 96 9.10502650 -0.001322229   1
## 4    3  0 0.06857102  0.416314640   1
## 5    4 22 3.20815215  0.114151375   1
## 6    5 86 8.95670223  0.434859137   1

We can also run multiple simulations simultaneously:

x <- simulate(ricker,nsim=10)
class(x)
## [1] "list"
sapply(x,class)
##  [1] "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp"
x <- simulate(ricker,nsim=10,as.data.frame=TRUE)
head(x)
##   time   y            N          e sim
## 1    0  61 7.000000e+00  0.0000000   1
## 2    1   3 2.400265e-01 -0.1729162   1
## 3    2  53 5.293097e+00 -0.4665640   1
## 4    3  17 1.443727e+00  0.1939209   1
## 5    4 155 1.868361e+01  0.2041458   1
## 6    5   0 8.847808e-06  0.3206253   1
str(x)
## 'data.frame':    510 obs. of  5 variables:
##  $ time: num  0 1 2 3 4 5 6 7 8 9 ...
##  $ y   : num  61 3 53 17 155 0 0 0 13 95 ...
##  $ N   : num  7 0.24 5.29 1.44 18.68 ...
##  $ e   : num  0 -0.173 -0.467 0.194 0.204 ...
##  $ sim : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...

Also,

x <- simulate(ricker,nsim=9,as.data.frame=TRUE,include.data=TRUE)
ggplot(data=x,aes(x=time,y=y,group=sim,color=(sim=="data")))+
  geom_line()+guides(color=FALSE)+
  facet_wrap(~sim,ncol=2)

We can compute a trajectory of the deterministic skeleton

y <- trajectory(ricker)
dim(y)
## [1]  2  1 51
dimnames(y)
## $variable
## [1] "N" "e"
## 
## $rep
## NULL
## 
## $time
## NULL
plot(time(ricker),y["N",1,],type="l")

Notice that ricker has parameters associated with it:

coef(ricker)
##        r    sigma      phi      N.0      e.0 
## 44.70118  0.30000 10.00000  7.00000  0.00000

These are the parameters at which the simulations and deterministic trajectory computations above were done. We can run these at different parameters:

theta <- coef(ricker)
theta[c("r","N.0")] <- c(5,3)
y <- trajectory(ricker,params=theta)
plot(time(ricker),y["N",1,],type="l")

x <- simulate(ricker,params=theta)
plot(x,var="y")

We can also change the parameters stored inside of ricker:

coef(ricker,c("r","N.0","sigma")) <- c(39,0.5,1)
coef(ricker)
##     r sigma   phi   N.0   e.0 
##  39.0   1.0  10.0   0.5   0.0
plot(simulate(ricker),var="y")

In all of the above, it’s possible to work with more than one set of parameters at a time. For example:

p <- parmat(coef(ricker),500)
dim(p); dimnames(p)
## [1]   5 500
## $variable
## [1] "r"     "sigma" "phi"   "N.0"   "e.0"  
## 
## $rep
## NULL
p["r",] <- seq(from=2,to=40,length=500)
y <- trajectory(ricker,params=p,times=200:1000)
matplot(p["r",],y["N",,],pch=".",col='black',xlab='r',ylab='N',log='x')

More information on manipulating and extracting information from pomp objects can be viewed in the help pages (methods?pomp).

There are a number of other examples included with the package. Do pompExamples() to see a list of these.

Inference algorithms in pomp

pomp provides a wide range of inference algorithms. We’ll learn about these in detail soon, but for now, let’s just look at some of their general features.

The pfilter function runs a simple particle filter. It can be used to evaluate the likelihood at a particular set of parameters. One uses the Np argument to specify the number of particles to use:

pf <- pfilter(ricker,Np=1000)
class(pf)
## [1] "pfilterd.pomp"
## attr(,"package")
## [1] "pomp"
plot(pf)

logLik(pf)
## [1] -200.4949

Note that pfilter returns an object of class pfilterd.pomp. This is the general rule: inference algorithms return objects that are pomp objects with additional information. The package provides tools to extract this information. We can run the particle filter again by doing

pf <- pfilter(pf)
logLik(pf)
## [1] -200.7841

which has the result of running the same computation again. Note that, because the particle filter is a Monte Carlo algorithm, we get a slightly different estimate of the log likelihood.

Note that, by default, running pfilter on a pfilterd.pomp object causes the computation to be re-run with the same parameters as before. Any additional arguments we add override these defaults. This is the general rule in pomp. For example,

pf <- pfilter(pf,Np=100)
logLik(pf)
## [1] -203.149

Here, the particle filtering has been performed with only 100 particles.

Building a custom pomp object

A real pomp data analysis begins with constructing one or more pomp objects to hold the data and the model or models under consideration. We’ll illustrate this process a dataset of Parus major abundance in Wytham Wood, near Oxford (McCleery and Perrins 1991).

Download and plot the data:

loc <- url("http://kinglab.eeb.lsa.umich.edu/SBIED/data/parus.csv")
dat <- read.csv(loc)
head(dat)
##   year pop
## 1 1960 148
## 2 1961 258
## 3 1962 185
## 4 1963 170
## 5 1964 267
## 6 1965 239
plot(pop~year,data=dat,type='o')

Let’s suppose that we want to fit the stochastic Ricker model discussed above to these data.

The call to construct a pomp object is, naturally enough, pomp. Documentation on this function can be had by doing ?pomp. Do class?pomp to get documentation on the pomp class. Learn about the various things you can do once you have a pomp object by doing methods?pomp and following the links therein. Read an overview of the package as a whole with links to its main features by doing package?pomp. A complete index of the functions in pomp is returned by the command library(help=pomp). Finally, the home page for the pomp project is ; there you have access to the complete source code, manuals, mailing lists, etc.

require(pomp)
parus <- pomp(dat,times="year",t0=1959)

The times argument specifies that the column labelled “year” gives the measurement times; t0 is the “zero-time”, the time at which the state process will be initialized. We’ve set it to one year prior to the beginning of the data. Plot it:

plot(parus)

Adding in the deterministic skeleton

We can add the Ricker model deterministic skeleton to the parus pomp object. Since the Ricker model is a discrete-time model, its skeleton is a map that takes \(N_t\) to \(N_{t+1}\) according to the Ricker model equation \[N_{t+1} = r\,N_{t}\,\exp(-N_{t}).\] We provide this to pomp in the form of a Csnippet, a little snippet of C code that performs the computation.

skel <- Csnippet("DN = r*N*exp(-N);")

We then add this to the pomp object:

parus <- pomp(parus,skeleton=skel,skeleton.type='map',
              paramnames=c("r"),statenames=c("N"))

Note that we have to inform pomp as to which of the variables we’ve referred to in skel is a state variable (statenames) and which is a parameter (paramnames). The skeleton.type argument tells pomp that the skeleton is a discrete-time dynamical system (a map) rather than a continuous-time system (a vectorfield).

With just the skeleton defined, we are in a position to compute the trajectories of the deterministic skeleton at any point in parameter space. For example,

traj <- trajectory(parus,params=c(N.0=1,r=12),as.data.frame=TRUE)
ggplot(data=traj,aes(x=time,y=N))+geom_line()

A note on terminology

If we know the state, \(x(t_0)\), of the system at time \(t_0\), it makes sense to speak about the entire trajectory of the system for all \(t>t_0\). This is true whether we are thinking of the system as deterministic or stochastic. Of course, in the former case, the trajectory is uniquely determined by \(x(t_0)\), while in the stochastic case, only the probability distribution of \(x(t)\), \(t>t_0\) is determined. In pomp, to avoid confusion, we use the term “trajectory” exclusively to refer to trajectories of a deterministic process. Thus, the trajectory command iterates or integrates the deterministic skeleton forward in time, returning the unique trajectory determined by the specified parameters. When we want to speak about sample paths of a stochastic process, we use the term simulation. Accordingly, the simulate command always returns individual sample paths from the POMP. In particular, we avoid “simulating a set of differential equations”, preferring instead to speak of “integrating” the equations, or “computing trajectories”.

Adding in the process model simulator

We can add the stochastic Ricker model to parus by writing a Csnippet that simulates one realization of the stochastic process, from an arbitary time \(t\) to \(t+1\), given arbitrary states and parameters. The following does this.

stochStep <- Csnippet("
  e = rnorm(0,sigma);
  N = r*N*exp(-N+e);
")
pomp(parus,rprocess=discrete.time.sim(step.fun=stochStep,delta.t=1),
     paramnames=c("r","sigma"),statenames=c("N","e")) -> parus

Note that in the above, we use the exp and rnorm functions from the R API. In general any C function provided by R is available to you. pomp also provides a number of C functions that are documented in the header file, pomp.h, that is installed with the package. See the Csnippet documentation (?Csnippet) to read more about how to write them. Note too that we use discrete.time.sim here because the model is a stochastic map. We specify that the time step of the discrete-time process is delta.t, here, 1 yr.

At this point, we have what we need to simulate the stochastic Ricker model.

sim <- simulate(parus,params=c(N.0=1,e.0=0,r=12,sigma=0.5),
                as.data.frame=TRUE,states=TRUE)
plot(N~time,data=sim,type='o')
lines(N~time,data=traj,type='l',col='red')

Adding in the measurement model and parameters

We complete the specification of the POMP by specifying the measurement model. To obtain the Poisson measurement model described above, we write two Csnippets. The first simulates:

rmeas <- Csnippet("pop = rpois(phi*N);")

and the second computes the likelihood of observing pop birds given a true density of N:

dmeas <- Csnippet("lik = dpois(pop,phi*N,give_log);")

[Note the give_log argument. When this code is evaluated, give_log will be set to 1 if the log likelihood is desired, and 0 else.] We add these into the pomp object:

pomp(parus,rmeasure=rmeas,dmeasure=dmeas,statenames=c("N"),paramnames=c("phi")) -> parus

Now we can simulate the whole POMP. First, let’s add some parameters to the pomp object:

coef(parus) <- c(N.0=1,e.0=0,r=20,sigma=0.1,phi=200)
sims <- simulate(parus,nsim=3,as.data.frame=TRUE,include.data=TRUE)
ggplot(data=sims,mapping=aes(x=time,y=pop))+geom_line()+
  facet_wrap(~sim)

Exercises

Ricker model parameters

Fiddle with the parameters to try and make the simulations look more like the data. This will help you build some intuition for what the various parameters do.

Reformulating the Ricker model

Reparameterize the Ricker model so that the scaling of \(N\) is explicit: \[N_{t+1} = r\,N_{t}\,\exp\left(-\frac{N_{t}}{K}\right).\]

Modify the pomp object we created above to reflect this reparameterization.

Modify the measurement model so that \[\mathrm{pop}_t \sim {\mathrm{Negbin}\left(\phi\,N_t,k\right)},\] i.e., \(\mathrm{pop}_t\) is negative-binomially distributed with mean \(\phi\,N_t\) and clumping parameter \(k\). See ?NegBinomial for documentation on the negative binomial distribution and the R Extensions Manual section on distribution functions for information on how to access these in C.

Beverton-Holt

Construct a pomp object for the Parus major data and the stochastic Beverton-Holt model \[N_{t+1} = \frac{a\,N_t}{1+b\,N_t}\,\varepsilon_t,\] where \(a\) and \(b\) are parameters and \[\varepsilon_t \sim {\mathrm{Lognormal}\left(-\tfrac{1}{2}\sigma^2,\sigma\right)}.\] Assume the same measurement model as before.

References

Bjørnstad, O. N., and B. T. Grenfell. 2001. Noisy clockwork: Time series analysis of population fluctuations in animals. Science 293:638–643.

He, D., E. L. Ionides, and A. A. King. 2010. Plug-and-play inference for disease dynamics: Measles in large and small populations as a case study. Journal of The Royal Society Interface 7:271–283.

King, A. A., M. Domenech de Cellès, F. M. G. Magpantay, and P. Rohani. 2015. Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to ebola. Proceedings of the Royal Society of London. Series B 282:20150347.

King, A. A., E. L. Ionides, M. Pascual, and M. J. Bouma. 2008. Inapparent infections and cholera dynamics. Nature 454:877–880.

Laneri, K., A. Bhadra, E. L. Ionides, M. Bouma, R. C. Dhiman, R. S. Yadav, and M. Pascual. 2010. Forcing versus feedback: Epidemic malaria and monsoon rains in northwest India. PLoS Computational Biology 6:e1000898.

McCleery, R. H., and C. M. Perrins. 1991. Effects of predation on the numbers of great tits Parus major. Pages 129–147 in Bird population studies. relevence to conservation and management. Oxford University Press, Oxford.

Romero-Severson, E. O., E. Volz, J. S. Koopman, T. Leitner, and E. L. Ionides. 2015. Dynamic variation in sexual contact rates in a cohort of HIV-negative gay men. American Journal of Epidemiology:kwv044.